Ch 16 exercises

The aim of these exercises is to continue making sure you’re comfortable with handling multivariate data. In this chapter’s, you’ll focus on analyses based on dissimilarities/distances, including fitting linear models to these kinds of response variables.

For each of the examples below, you should follow the sequence we’ve used previously, as far as it’s sensible:

1. What is the biological question?

2. Is the predictor continuous or categorical?

3. Write out the linear model corresponding to this question.

4. What distribution do you expect the response variable to follow?

5. What are the assumptions behind the statistical model you’ll fit?

   1. Are those assumptions satisfied?

6. Fit the model

   1. How will you assess whether the model fits well?

   2. Can you detect an effect of the predictors?

   3. How do you measure the effect?

7. What do you conclude (including any cautions)

---

A

Dixon et al. (2018) focused on assemblages of reptiles in forests and woodlands of southeastern Australia, with a particular interest in relationships between these assemblages and the fire history of the landscape. They identified 81 sites that varied in time since fire, from 6 months to >96 y. Rather than a continuum, three categories of time since fire were used, 0.5-2y, 6-12y, and >96y. Sites were also classified according to habitat.

Reptile assemblages were sampled with a range of methods, including visual surveys and camera traps, to give counts of 20 reptiles, whose abundance ranged from 0-1 to 0-126, depending on species.

Data are available from dryad, as Rep_abund.csv. You’ll want to focus on the columns with reptile numbers, and the one with the fire history (tsf). For convenience, we’ve extracted those data as dixonbiota.csv

Start by having a quick look at the data.

Rows: 79 Columns: 22── Column specification ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr  (2): site, tsf
dbl (20): adup, aplat, amur, amac, aram, dcor, ecun, esax, eul, htal, ldel, lgui, lwhi, ppor, pent, pspe, ptex, rdie, tnig, vros
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

The file is pretty straightforward, with tsf being the fire category, and columns 3-22 each representing one reptile taxon.

Outline how you’d assess whether the reptile assemblage (the combination of species present and their abundances) is related to fire history, using a dissimilarity-based approach.

You should think about how you’ll deal with the data:

• Will you need a transformation, standardization, etc.?

• What are the consequences of any decisions you make for the interpretation of your analysis?

• What measure(s) of dissimilarity are appropriate?

Run the analysis and provide your interpretation

MDS on reptile abundances

Start with dissimilarities

Students should have some discussion of raw vs transformed (e.g. 4th root), standardization.

Fourth root is common. Using it downweights the influence of abundant species, while raw allows those species to have a strong influence. There are questions about which one of these reflects ecological function better

Data are all abundances

Visualize results

Need to decide whether to use 2 or 3 dimensions

Run 0 stress 0.184364 
Run 1 stress 0.1954244 
Run 2 stress 0.1835248 
... New best solution
... Procrustes: rmse 0.02400721  max resid 0.1627118 
Run 3 stress 0.1875543 
Run 4 stress 0.1873871 
Run 5 stress 0.1835255 
... Procrustes: rmse 0.0002431954  max resid 0.001211121 
... Similar to previous best
Run 6 stress 0.1835255 
... Procrustes: rmse 0.0004322375  max resid 0.002396493 
... Similar to previous best
Run 7 stress 0.187394 
Run 8 stress 0.1839905 
... Procrustes: rmse 0.01513281  max resid 0.1010331 
Run 9 stress 0.1835253 
... Procrustes: rmse 0.0003986974  max resid 0.002827447 
... Similar to previous best
Run 10 stress 0.1875557 
Run 11 stress 0.1835394 
... Procrustes: rmse 0.00239692  max resid 0.01890764 
Run 12 stress 0.1875561 
Run 13 stress 0.1962511 
Run 14 stress 0.1845587 
Run 15 stress 0.1835253 
... Procrustes: rmse 0.0007731379  max resid 0.004859597 
... Similar to previous best
Run 16 stress 0.1855483 
Run 17 stress 0.1901239 
Run 18 stress 0.1876739 
Run 19 stress 0.1874628 
Run 20 stress 0.1837457 
... Procrustes: rmse 0.009486057  max resid 0.06288526 
Run 21 stress 0.1835247 
... New best solution
... Procrustes: rmse 0.000444541  max resid 0.003456569 
... Similar to previous best
Run 22 stress 0.1835247 
... Procrustes: rmse 0.000423645  max resid 0.0032015 
... Similar to previous best
Run 23 stress 0.1893602 
Run 24 stress 0.1873867 
Run 25 stress 0.1853823 
Run 26 stress 0.1923026 
Run 27 stress 0.1839242 
... Procrustes: rmse 0.01343186  max resid 0.08937936 
Run 28 stress 0.1880216 
Run 29 stress 0.1875504 
Run 30 stress 0.1844798 
Run 31 stress 0.1835267 
... Procrustes: rmse 0.001099299  max resid 0.00825338 
... Similar to previous best
Run 32 stress 0.1880227 
Run 33 stress 0.1880842 
Run 34 stress 0.1874227 
Run 35 stress 0.1835265 
... Procrustes: rmse 0.000828435  max resid 0.006610793 
... Similar to previous best
Run 36 stress 0.1835258 
... Procrustes: rmse 0.0009260963  max resid 0.007366734 
... Similar to previous best
Run 37 stress 0.1875538 
Run 38 stress 0.1880843 
Run 39 stress 0.1844375 
Run 40 stress 0.1835246 
... New best solution
... Procrustes: rmse 5.771551e-05  max resid 0.0002761158 
... Similar to previous best
*** Best solution repeated 1 times

Call:
metaMDS(comm = df1s.bc, k = 3, try = 40, trymax = 80, autotransform = FALSE,      maxit = 200) 

global Multidimensional Scaling using monoMDS

Data:     df1s.bc 
Distance: bray 

Dimensions: 3 
Stress:     0.1835246 
Stress type 1, weak ties
Best solution was repeated 1 time in 40 tries
The best solution was from try 40 (random start)
Scaling: centring, PC rotation, halfchange scaling 
Species: scores missing

Run 0 stress 0.2620287 
Run 1 stress 0.2657228 
Run 2 stress 0.2613208 
... New best solution
... Procrustes: rmse 0.03050854  max resid 0.1997569 
Run 3 stress 0.2634799 
Run 4 stress 0.255339 
... New best solution
... Procrustes: rmse 0.08914296  max resid 0.2703881 
Run 5 stress 0.2585276 
Run 6 stress 0.2557278 
... Procrustes: rmse 0.008721101  max resid 0.04949448 
Run 7 stress 0.2557123 
... Procrustes: rmse 0.007310089  max resid 0.03878651 
Run 8 stress 0.2553666 
... Procrustes: rmse 0.0102181  max resid 0.07939793 
Run 9 stress 0.2700152 
Run 10 stress 0.2644649 
Run 11 stress 0.2626207 
Run 12 stress 0.2629297 
Run 13 stress 0.2650439 
Run 14 stress 0.268519 
Run 15 stress 0.2684536 
Run 16 stress 0.2622278 
Run 17 stress 0.2692256 
Run 18 stress 0.268778 
Run 19 stress 0.2643455 
Run 20 stress 0.2607326 
Run 21 stress 0.2685011 
Run 22 stress 0.2637274 
Run 23 stress 0.2652818 
Run 24 stress 0.2717004 
Run 25 stress 0.2597301 
Run 26 stress 0.258373 
Run 27 stress 0.2629686 
Run 28 stress 0.2736616 
Run 29 stress 0.2573166 
Run 30 stress 0.3044425 
Run 31 stress 0.2715289 
Run 32 stress 0.2585135 
Run 33 stress 0.2618396 
Run 34 stress 0.2557593 
... Procrustes: rmse 0.01141604  max resid 0.06725798 
Run 35 stress 0.2567309 
Run 36 stress 0.2621995 
Run 37 stress 0.2556597 
... Procrustes: rmse 0.007557198  max resid 0.03638444 
Run 38 stress 0.2661423 
Run 39 stress 0.2621501 
Run 40 stress 0.2626279 
Run 41 stress 0.2638528 
Run 42 stress 0.2666209 
Run 43 stress 0.2552761 
... New best solution
... Procrustes: rmse 0.003378444  max resid 0.0214193 
Run 44 stress 0.2643959 
Run 45 stress 0.2559259 
Run 46 stress 0.2702403 
Run 47 stress 0.2743521 
Run 48 stress 0.2639517 
Run 49 stress 0.2674239 
Run 50 stress 0.2660627 
Run 51 stress 0.260206 
Run 52 stress 0.259709 
Run 53 stress 0.2695475 
Run 54 stress 0.2652297 
Run 55 stress 0.2669659 
Run 56 stress 0.2695288 
Run 57 stress 0.2585219 
Run 58 stress 0.2553393 
... Procrustes: rmse 0.003381346  max resid 0.02138202 
Run 59 stress 0.2635542 
Run 60 stress 0.2603403 
Run 61 stress 0.2665924 
Run 62 stress 0.2673204 
Run 63 stress 0.2627007 
Run 64 stress 0.2691563 
Run 65 stress 0.2634604 
Run 66 stress 0.2557524 
... Procrustes: rmse 0.01015447  max resid 0.06110719 
Run 67 stress 0.2587965 
Run 68 stress 0.2600645 
Run 69 stress 0.2673046 
Run 70 stress 0.2636262 
Run 71 stress 0.2629212 
Run 72 stress 0.2569404 
Run 73 stress 0.2587754 
Run 74 stress 0.2628626 
Run 75 stress 0.2662464 
Run 76 stress 0.2693245 
Run 77 stress 0.2701223 
Run 78 stress 0.2646634 
Run 79 stress 0.2667217 
Run 80 stress 0.2666245 
*** Best solution was not repeated -- monoMDS stopping criteria:
    13: no. of iterations >= maxit
    67: stress ratio > sratmax

Call:
metaMDS(comm = df1s.bc, k = 2, try = 40, trymax = 80, autotransform = FALSE,      maxit = 200) 

global Multidimensional Scaling using monoMDS

Data:     df1s.bc 
Distance: bray 

Dimensions: 2 
Stress:     0.2552761 
Stress type 1, weak ties
Best solution was not repeated after 80 tries
The best solution was from try 43 (random start)
Scaling: centring, PC rotation, halfchange scaling 
Species: scores missing

Stress values suggest 3d plot is a more accurate representation

Need to plot dim 1 vs dim 2 and dim 1 vs dim 3 with tsf groups identified by symbols

Generate graphs

This is a variant of code used in worked example for Ch 16. It does use the graph settings defined in appearance.R, which is called earlier.

How would you fit a linear model to assess fire effects?

Do permanova including tsf

Permutation test for adonis under reduced model
Terms added sequentially (first to last)
Permutation: free
Number of permutations: 999

adonis2(formula = df1s.bc ~ tsf, data = df, permutations = 999)
         Df SumOfSqs      R2      F Pr(>F)    
tsf       2   2.6433 0.11039 4.7151  0.001 ***
Residual 76  21.3028 0.88961                  
Total    78  23.9460 1.00000                  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Permanova suggests strong signal from tsf

SIMPER analysis for tsf group differences

Contrast: long_short 

     average      sd   ratio     ava     avb cumsum      p    
lgui 0.13247 0.14056 0.94250 0.15918 0.26922  0.172 0.3247    
eul  0.09332 0.12167 0.76700 0.02379 0.16793  0.293 0.0430 *  
dcor 0.08201 0.11054 0.74190 0.06582 0.10672  0.399 0.1219    
amac 0.06242 0.05874 1.06250 0.07712 0.07023  0.480 0.9760    
tnig 0.05948 0.07551 0.78760 0.03804 0.07852  0.557 0.8492    
ldel 0.05168 0.04157 1.24340 0.08687 0.00000  0.624 0.0010 ***
pent 0.05090 0.04718 1.07890 0.10975 0.03392  0.690 0.7133    
htal 0.03859 0.05635 0.68490 0.04599 0.02950  0.740 0.4176    
lwhi 0.03208 0.03938 0.81470 0.04549 0.01261  0.782 0.0160 *  
pspe 0.02827 0.03740 0.75580 0.04834 0.00680  0.819 0.0539 .  
amur 0.02647 0.05251 0.50400 0.04512 0.00000  0.853 0.9171    
rdie 0.02561 0.09803 0.26120 0.00512 0.04275  0.886 0.5365    
ecun 0.02190 0.04456 0.49150 0.03603 0.00000  0.915 0.0420 *  
esax 0.01860 0.04544 0.40940 0.01674 0.01687  0.939 0.1259    
ptex 0.01585 0.04231 0.37450 0.01170 0.01687  0.959 0.2318    
aram 0.01495 0.03676 0.40680 0.01517 0.01261  0.979 0.5894    
vros 0.01279 0.03799 0.33650 0.02216 0.00000  0.995 0.8422    
ppor 0.00363 0.01891 0.19220 0.00620 0.00000  1.000 0.7063    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Contrast: long_medium 

     average      sd   ratio     ava     avb cumsum      p    
lgui 0.10221 0.10333 0.98920 0.15918 0.21705  0.133 0.9720    
amac 0.10139 0.11774 0.86110 0.07712 0.16122  0.264 0.2078    
tnig 0.07029 0.10274 0.68410 0.03804 0.10788  0.356 0.6713    
pent 0.06601 0.06862 0.96190 0.10975 0.08112  0.442 0.0260 *  
amur 0.05952 0.11624 0.51210 0.04512 0.07283  0.519 0.0649 .  
dcor 0.05521 0.09128 0.60480 0.06582 0.03801  0.591 0.8032    
ldel 0.04942 0.03935 1.25580 0.08687 0.00185  0.655 0.0010 ***
eul  0.04752 0.07141 0.66540 0.02379 0.07484  0.716 0.9970    
htal 0.03907 0.06803 0.57430 0.04599 0.02956  0.767 0.3227    
vros 0.03267 0.09018 0.36230 0.02216 0.03599  0.809 0.2198    
lwhi 0.02992 0.03771 0.79360 0.04549 0.00840  0.848 0.0010 ***
pspe 0.02952 0.03909 0.75530 0.04834 0.00806  0.887 0.0010 ***
ecun 0.02116 0.04260 0.49670 0.03603 0.00000  0.914 0.0020 ** 
rdie 0.01844 0.06391 0.28850 0.00512 0.02970  0.938 0.8292    
aram 0.01711 0.04409 0.38800 0.01517 0.01739  0.960 0.4705    
esax 0.01138 0.02737 0.41600 0.01674 0.00420  0.975 0.6653    
ptex 0.00983 0.02516 0.39070 0.01170 0.00658  0.988 0.7073    
ppor 0.00929 0.03741 0.24840 0.00620 0.01107  1.000 0.3616    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Contrast: short_medium 

     average      sd   ratio     ava     avb cumsum      p   
lgui 0.15998 0.16074 0.99530 0.26922 0.21705  0.205 0.0140 * 
amac 0.10618 0.14346 0.74010 0.07023 0.16122  0.340 0.1738   
eul  0.10520 0.12391 0.84900 0.16793 0.07484  0.475 0.0070 **
tnig 0.08903 0.11469 0.77630 0.07852 0.10788  0.589 0.0799 . 
dcor 0.07378 0.14450 0.51060 0.10672 0.03801  0.683 0.2388   
pent 0.04966 0.08217 0.60440 0.03392 0.08112  0.747 0.7692   
amur 0.04038 0.12310 0.32800 0.00000 0.07283  0.798 0.7073   
rdie 0.03751 0.11223 0.33420 0.04275 0.02970  0.847 0.1508   
htal 0.03122 0.07878 0.39630 0.02950 0.02956  0.886 0.6993   
vros 0.02238 0.09233 0.24240 0.00000 0.03599  0.915 0.6074   
aram 0.01559 0.04775 0.32650 0.01261 0.01739  0.935 0.5524   
ptex 0.01278 0.04462 0.28650 0.01687 0.00658  0.951 0.4316   
esax 0.01161 0.04217 0.27530 0.01687 0.00420  0.966 0.5704   
lwhi 0.01095 0.03778 0.28970 0.01261 0.00840  0.980 0.9900   
pspe 0.00840 0.02396 0.35050 0.00680 0.00806  0.991 0.9980   
ppor 0.00608 0.03495 0.17380 0.00000 0.01107  0.999 0.5385   
ldel 0.00101 0.00583 0.17380 0.00000 0.00185  1.000 1.0000   
ecun 0.00000 0.00000     NaN 0.00000 0.00000  1.000 0.9990   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Permutation: free
Number of permutations: 1000

Suppose you think of a reptile assemblage as simply the species that are present, regardless of how common they are. What would you conclude about the infuence of time since fire now?

Repeat the preceding analysis after creating a new data file that is presence-absence only.

Run 0 stress 0.1404128 
Run 1 stress 0.140413 
... Procrustes: rmse 0.000158617  max resid 0.0009441077 
... Similar to previous best
Run 2 stress 0.1489204 
Run 3 stress 0.1448005 
Run 4 stress 0.1486182 
Run 5 stress 0.1453965 
Run 6 stress 0.1418091 
Run 7 stress 0.1415816 
Run 8 stress 0.1481748 
Run 9 stress 0.1418881 
Run 10 stress 0.1404115 
... New best solution
... Procrustes: rmse 0.0004234926  max resid 0.002436781 
... Similar to previous best
Run 11 stress 0.1488108 
Run 12 stress 0.1482426 
Run 13 stress 0.1419245 
Run 14 stress 0.1590438 
Run 15 stress 0.143286 
Run 16 stress 0.1486203 
Run 17 stress 0.1405587 
... Procrustes: rmse 0.01449226  max resid 0.08640004 
Run 18 stress 0.1486316 
Run 19 stress 0.1405592 
... Procrustes: rmse 0.01479786  max resid 0.08730948 
Run 20 stress 0.1405595 
... Procrustes: rmse 0.01483792  max resid 0.08736059 
Run 21 stress 0.1404114 
... New best solution
... Procrustes: rmse 0.002731297  max resid 0.02192069 
Run 22 stress 0.1496581 
Run 23 stress 0.1416846 
Run 24 stress 0.1405862 
... Procrustes: rmse 0.01597208  max resid 0.09180616 
Run 25 stress 0.1461628 
Run 26 stress 0.1404123 
... Procrustes: rmse 0.000245968  max resid 0.001978068 
... Similar to previous best
Run 27 stress 0.1405596 
... Procrustes: rmse 0.01457381  max resid 0.08827165 
Run 28 stress 0.1404116 
... Procrustes: rmse 0.0002275687  max resid 0.001455665 
... Similar to previous best
Run 29 stress 0.1478219 
Run 30 stress 0.1438487 
Run 31 stress 0.1404119 
... Procrustes: rmse 0.0001425049  max resid 0.001044504 
... Similar to previous best
Run 32 stress 0.1478221 
Run 33 stress 0.141641 
Run 34 stress 0.1414341 
Run 35 stress 0.1422857 
Run 36 stress 0.1487658 
Run 37 stress 0.1486133 
Run 38 stress 0.1502159 
Run 39 stress 0.1486756 
Run 40 stress 0.1404122 
... Procrustes: rmse 0.0002283066  max resid 0.001565205 
... Similar to previous best
*** Best solution repeated 4 times

Call:
metaMDS(comm = df1.jac, k = 3, try = 40, trymax = 80, autotransform = FALSE,      maxit = 200) 

global Multidimensional Scaling using monoMDS

Data:     df1.jac 
Distance: binary jaccard 

Dimensions: 3 
Stress:     0.1404114 
Stress type 1, weak ties
Best solution was repeated 4 times in 40 tries
The best solution was from try 21 (random start)
Scaling: centring, PC rotation, halfchange scaling 
Species: scores missing

Run 0 stress 0.2009509 
Run 1 stress 0.2019149 
Run 2 stress 0.2011657 
... Procrustes: rmse 0.006416695  max resid 0.0398645 
Run 3 stress 0.202576 
Run 4 stress 0.2284835 
Run 5 stress 0.2197434 
Run 6 stress 0.2011646 
... Procrustes: rmse 0.006850207  max resid 0.03839815 
Run 7 stress 0.2204112 
Run 8 stress 0.215058 
Run 9 stress 0.203002 
Run 10 stress 0.2143484 
Run 11 stress 0.2104918 
Run 12 stress 0.2180894 
Run 13 stress 0.2279185 
Run 14 stress 0.2071546 
Run 15 stress 0.2160939 
Run 16 stress 0.2164477 
Run 17 stress 0.2108417 
Run 18 stress 0.2097709 
Run 19 stress 0.2079378 
Run 20 stress 0.2138375 
Run 21 stress 0.2120034 
Run 22 stress 0.2115225 
Run 23 stress 0.215718 
Run 24 stress 0.2057211 
Run 25 stress 0.2030018 
Run 26 stress 0.213385 
Run 27 stress 0.2285741 
Run 28 stress 0.2289543 
Run 29 stress 0.2178744 
Run 30 stress 0.208991 
Run 31 stress 0.2135402 
Run 32 stress 0.20361 
Run 33 stress 0.2156126 
Run 34 stress 0.2146268 
Run 35 stress 0.2167606 
Run 36 stress 0.2274825 
Run 37 stress 0.202578 
Run 38 stress 0.2036086 
Run 39 stress 0.2171294 
Run 40 stress 0.2366921 
Run 41 stress 0.2011785 
... Procrustes: rmse 0.006149074  max resid 0.03876166 
Run 42 stress 0.218512 
Run 43 stress 0.2020987 
Run 44 stress 0.2068717 
Run 45 stress 0.2041772 
Run 46 stress 0.203832 
Run 47 stress 0.2272148 
Run 48 stress 0.2035363 
Run 49 stress 0.2074059 
Run 50 stress 0.2174067 
Run 51 stress 0.2098787 
Run 52 stress 0.2074132 
Run 53 stress 0.2011648 
... Procrustes: rmse 0.006869384  max resid 0.03837032 
Run 54 stress 0.2133796 
Run 55 stress 0.2139331 
Run 56 stress 0.2112664 
Run 57 stress 0.226952 
Run 58 stress 0.2011782 
... Procrustes: rmse 0.006123098  max resid 0.03888408 
Run 59 stress 0.2119751 
Run 60 stress 0.2133246 
Run 61 stress 0.2134668 
Run 62 stress 0.2045553 
Run 63 stress 0.2084428 
Run 64 stress 0.2019056 
Run 65 stress 0.2041934 
Run 66 stress 0.2144476 
Run 67 stress 0.2068721 
Run 68 stress 0.2128513 
Run 69 stress 0.2069873 
Run 70 stress 0.2116038 
Run 71 stress 0.2187148 
Run 72 stress 0.2030068 
Run 73 stress 0.2078151 
Run 74 stress 0.2011768 
... Procrustes: rmse 0.005962949  max resid 0.03913871 
Run 75 stress 0.2056871 
Run 76 stress 0.2011778 
... Procrustes: rmse 0.006066117  max resid 0.03889517 
Run 77 stress 0.2177676 
Run 78 stress 0.2127882 
Run 79 stress 0.2111792 
Run 80 stress 0.2014982 
*** Best solution was not repeated -- monoMDS stopping criteria:
    12: no. of iterations >= maxit
    68: stress ratio > sratmax

Call:
metaMDS(comm = df1.jac, k = 2, try = 40, trymax = 80, autotransform = FALSE,      maxit = 200) 

global Multidimensional Scaling using monoMDS

Data:     df1.jac 
Distance: binary jaccard 

Dimensions: 2 
Stress:     0.2009509 
Stress type 1, weak ties
Best solution was not repeated after 80 tries
The best solution was from try 0 (metric scaling or null solution)
Scaling: centring, PC rotation, halfchange scaling 
Species: scores missing

Marginal decision here - stress just over 0.2 for k=2, so might be OK. Stress good for k=3, but 2-d MDS usually better for reporting or showing to audience

Unburned sites separate very clearly on MDS1. Pattern clearer than with abundances included

Do permanova including tsf

Permutation test for adonis under reduced model
Terms added sequentially (first to last)
Permutation: free
Number of permutations: 999

adonis2(formula = df1.jac ~ tsf, data = df, permutations = 999)
         Df SumOfSqs      R2      F Pr(>F)    
tsf       2   3.7853 0.20678 9.9057  0.001 ***
Residual 76  14.5210 0.79322                  
Total    78  18.3062 1.00000                  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Clear fire signal

Extension activity

Dixon and her colleagues also recorded a range of habitat variables, which we’ve collected for you in dixonenv.csv

Rows: 79 Columns: 14── Column specification ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr  (4): site, tsf, veg, aspect
dbl (10): lcwd, shrcov, grcov, litcov, litdep, rocks, elev, warm, cold, twi
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

They recorded Coarse Woody Debris, which was long-transformed, litter cover (litcov), % cover of groundcover (grcov), and % cover of shrubs (shrcov), and rock cover. In the original paper, Dixon et al. were comfortable that these predictors weren’t correlated.

How is reptile assemblage related to time since fire and these five habitat variables

Check continuous predictors

Analysis of Variance Table

Response: Distances
          Df  Sum Sq  Mean Sq F value    Pr(>F)    
Groups     2 0.48563 0.242815  32.165 7.573e-11 ***
Residuals 76 0.57374 0.007549                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

do permanova including tsf plus continuous predictors

Permutation test for adonis under reduced model
Terms added sequentially (first to last)
Permutation: free
Number of permutations: 999

adonis2(formula = df1s.bc ~ tsf + lcwd + litcov + grcov + shrcov + lrocks, data = dixonenv, permutations = 999)
         Df SumOfSqs      R2      F Pr(>F)    
tsf       2   2.6433 0.11039 4.8306  0.001 ***
lcwd      1   0.4865 0.02032 1.7783  0.055 .  
litcov    1   0.4037 0.01686 1.4757  0.128    
grcov     1   0.3281 0.01370 1.1991  0.284    
shrcov    1   0.3291 0.01374 1.2028  0.253    
lrocks    1   0.3300 0.01378 1.2062  0.274    
Residual 71  19.4253 0.81121                  
Total    78  23.9460 1.00000                  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Permutation test for adonis under reduced model
Terms added sequentially (first to last)
Permutation: free
Number of permutations: 999

adonis2(formula = df1.jac ~ tsf + lcwd + litcov + grcov + shrcov + lrocks, data = dixonenv, permutations = 999)
         Df SumOfSqs      R2       F Pr(>F)    
tsf       2   3.7853 0.20678 10.4082  0.001 ***
lcwd      1   0.3802 0.02077  2.0908  0.032 *  
litcov    1   0.3145 0.01718  1.7296  0.072 .  
grcov     1   0.4097 0.02238  2.2530  0.015 *  
shrcov    1   0.2748 0.01501  1.5115  0.124    
lrocks    1   0.2309 0.01262  1.2700  0.238    
Residual 71  12.9108 0.70527                   
Total    78  18.3062 1.00000                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Main effect is tsf; only some question about CWD when we look at species composition and abundance, while for presence-absence data, there are also effects of CWD and groundcover

mvabund

Error in model.frame.default(formula = dfmv ~ dixonenv$tsf + dixonenv$lcwd +  : 
  variable lengths differ (found for 'dixonenv$tsf')

B Simple MDS & Permanova

Let’s return to the Hutto and Barrett (2021) example from the Chapter 15 exercises. There, we used similarity-based analyses to explore the relationship between frog assemblages in ponds and the degree of urbanization surrounding. This seems a question that could just as easily be examined using distances or dissimilarities.

Return to the data and assess whether frog assemblages (using the standardized abundance scale or presence-absence) differ with urbanization.

Did your conclusions differ from those obtained using RDA?

Rows: 50 Columns: 14── Column specification ──────────────────────────────────────────────────────────────────────
Delimiter: ","
chr  (1): type
dbl (13): Site, ACRE, BAME, BFOW, GCAR, HCIN, HVER, LCAT, LCLA, LPAL, LSPH, PCRU, PFER
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

The explanation of the variables is in the previous chapter’s exercises.

Run 0 stress 0.196708 
Run 1 stress 0.1877627 
... New best solution
... Procrustes: rmse 0.06210972  max resid 0.2289889 
Run 2 stress 0.1965252 
Run 3 stress 0.2177645 
Run 4 stress 0.1983645 
Run 5 stress 0.2034451 
Run 6 stress 0.2052112 
Run 7 stress 0.1877626 
... New best solution
... Procrustes: rmse 0.0001156879  max resid 0.0003505044 
... Similar to previous best
Run 8 stress 0.185584 
... New best solution
... Procrustes: rmse 0.07743106  max resid 0.5021775 
Run 9 stress 0.1992463 
Run 10 stress 0.1974204 
Run 11 stress 0.1946875 
Run 12 stress 0.1862779 
Run 13 stress 0.1863292 
Run 14 stress 0.2005585 
Run 15 stress 0.1920873 
Run 16 stress 0.1875597 
Run 17 stress 0.1964527 
Run 18 stress 0.1920872 
Run 19 stress 0.1977907 
Run 20 stress 0.2167665 
Run 21 stress 0.1875595 
Run 22 stress 0.1929107 
Run 23 stress 0.1975661 
Run 24 stress 0.186278 
Run 25 stress 0.1875598 
Run 26 stress 0.1960212 
Run 27 stress 0.1991282 
Run 28 stress 0.1975497 
Run 29 stress 0.185397 
... New best solution
... Procrustes: rmse 0.02114207  max resid 0.1079972 
Run 30 stress 0.1969418 
Run 31 stress 0.186278 
Run 32 stress 0.1862238 
Run 33 stress 0.1877806 
Run 34 stress 0.1919384 
Run 35 stress 0.1967647 
Run 36 stress 0.2073798 
Run 37 stress 0.1888207 
Run 38 stress 0.2054248 
Run 39 stress 0.1874611 
Run 40 stress 0.1862237 
Run 41 stress 0.186278 
Run 42 stress 0.1888251 
Run 43 stress 0.1937375 
Run 44 stress 0.1862237 
Run 45 stress 0.1963067 
Run 46 stress 0.1862237 
Run 47 stress 0.1854373 
... Procrustes: rmse 0.009296219  max resid 0.0491063 
Run 48 stress 0.1920872 
Run 49 stress 0.1875595 
Run 50 stress 0.1853641 
... New best solution
... Procrustes: rmse 0.001953553  max resid 0.00909663 
... Similar to previous best
*** Best solution repeated 1 times

Call:
metaMDS(comm = df1.bc, k = 2, try = 40, trymax = 80, autotransform = FALSE,      maxit = 200) 

global Multidimensional Scaling using monoMDS

Data:     df1.bc 
Distance: bray 

Dimensions: 2 
Stress:     0.1853641 
Stress type 1, weak ties
Best solution was repeated 1 time in 50 tries
The best solution was from try 50 (random start)
Scaling: centring, PC rotation, halfchange scaling 
Species: scores missing

Run 0 stress 0.129187 
Run 1 stress 0.1291859 
... New best solution
... Procrustes: rmse 0.0008388266  max resid 0.00496375 
... Similar to previous best
Run 2 stress 0.1402078 
Run 3 stress 0.1308799 
Run 4 stress 0.1289613 
... New best solution
... Procrustes: rmse 0.01088561  max resid 0.05621761 
Run 5 stress 0.129186 
... Procrustes: rmse 0.01089164  max resid 0.05582815 
Run 6 stress 0.1289616 
... Procrustes: rmse 0.0002674815  max resid 0.001155959 
... Similar to previous best
Run 7 stress 0.1380881 
Run 8 stress 0.1425061 
Run 9 stress 0.1395557 
Run 10 stress 0.1291859 
... Procrustes: rmse 0.01090439  max resid 0.05603889 
Run 11 stress 0.1291861 
... Procrustes: rmse 0.01047698  max resid 0.0544579 
Run 12 stress 0.1380864 
Run 13 stress 0.12921 
... Procrustes: rmse 0.01243608  max resid 0.06398953 
Run 14 stress 0.1291867 
... Procrustes: rmse 0.01108112  max resid 0.05658361 
Run 15 stress 0.1380875 
Run 16 stress 0.1289612 
... New best solution
... Procrustes: rmse 0.0001897893  max resid 0.001008599 
... Similar to previous best
Run 17 stress 0.1291861 
... Procrustes: rmse 0.01047209  max resid 0.05451001 
Run 18 stress 0.1384298 
Run 19 stress 0.1291858 
... Procrustes: rmse 0.01085358  max resid 0.05610837 
Run 20 stress 0.1289617 
... Procrustes: rmse 0.0005107488  max resid 0.002356854 
... Similar to previous best
Run 21 stress 0.129192 
... Procrustes: rmse 0.01049379  max resid 0.05482609 
Run 22 stress 0.1291861 
... Procrustes: rmse 0.01090769  max resid 0.05651704 
Run 23 stress 0.1291858 
... Procrustes: rmse 0.01084092  max resid 0.05603118 
Run 24 stress 0.1291603 
... Procrustes: rmse 0.009742854  max resid 0.05179386 
Run 25 stress 0.1291861 
... Procrustes: rmse 0.01048144  max resid 0.05456092 
Run 26 stress 0.1291857 
... Procrustes: rmse 0.01077451  max resid 0.05553841 
Run 27 stress 0.1380876 
Run 28 stress 0.1402635 
Run 29 stress 0.1384275 
Run 30 stress 0.1291864 
... Procrustes: rmse 0.01045087  max resid 0.05448513 
Run 31 stress 0.129186 
... Procrustes: rmse 0.01090378  max resid 0.0564044 
Run 32 stress 0.1291859 
... Procrustes: rmse 0.01083027  max resid 0.0561575 
Run 33 stress 0.1380913 
Run 34 stress 0.1291855 
... Procrustes: rmse 0.01069769  max resid 0.05541034 
Run 35 stress 0.1383999 
Run 36 stress 0.1289618 
... Procrustes: rmse 0.0005305334  max resid 0.002504772 
... Similar to previous best
Run 37 stress 0.1291859 
... Procrustes: rmse 0.01084244  max resid 0.05590108 
Run 38 stress 0.1289621 
... Procrustes: rmse 0.0006163776  max resid 0.002721091 
... Similar to previous best
Run 39 stress 0.1291858 
... Procrustes: rmse 0.01057936  max resid 0.05497957 
Run 40 stress 0.1291858 
... Procrustes: rmse 0.01050048  max resid 0.05461872 
*** Best solution repeated 4 times

Call:
metaMDS(comm = df1.bc, k = 3, try = 40, trymax = 80, autotransform = FALSE,      maxit = 200) 

global Multidimensional Scaling using monoMDS

Data:     df1.bc 
Distance: bray 

Dimensions: 3 
Stress:     0.1289612 
Stress type 1, weak ties
Best solution was repeated 4 times in 40 tries
The best solution was from try 16 (random start)
Scaling: centring, PC rotation, halfchange scaling 
Species: scores missing

2-d is acceptable

Permutation test for adonis under reduced model
Terms added sequentially (first to last)
Permutation: free
Number of permutations: 999

adonis2(formula = df1.bc ~ type, data = df, permutations = 999)
         Df SumOfSqs      R2      F Pr(>F)  
type      2    0.912 0.07281 1.7669  0.034 *
Residual 45   11.614 0.92719                
Total    47   12.526 1.00000                
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Permanova detects a difference; MDS plot suggests low urbanization may be different from the other two, though worth looking at dispersion as well, as spread of this type is less than the other two on plot

Repeat analysis with presence-absence

df1.jac <- vegdist(df1,binary=TRUE,method=‘jaccard’)

Run 0 stress 0.1551612 
Run 1 stress 0.1538244 
... New best solution
... Procrustes: rmse 0.01653307  max resid 0.06524188 
Run 2 stress 0.1808491 
Run 3 stress 0.1646706 
Run 4 stress 0.1901198 
Run 5 stress 0.1908161 
Run 6 stress 0.1607835 
Run 7 stress 0.1961593 
Run 8 stress 0.163888 
Run 9 stress 0.158713 
Run 10 stress 0.1551611 
Run 11 stress 0.1728782 
Run 12 stress 0.1627091 
Run 13 stress 0.1602264 
Run 14 stress 0.1538803 
... Procrustes: rmse 0.0179707  max resid 0.07573361 
Run 15 stress 0.1609335 
Run 16 stress 0.1646707 
Run 17 stress 0.1945193 
Run 18 stress 0.1837141 
Run 19 stress 0.1669303 
Run 20 stress 0.1585798 
Run 21 stress 0.1630282 
Run 22 stress 0.15858 
Run 23 stress 0.1585799 
Run 24 stress 0.1734741 
Run 25 stress 0.1630923 
Run 26 stress 0.1585799 
Run 27 stress 0.1602583 
Run 28 stress 0.15858 
Run 29 stress 0.1752544 
Run 30 stress 0.1784224 
Run 31 stress 0.1609031 
Run 32 stress 0.1553998 
Run 33 stress 0.1707427 
Run 34 stress 0.1604495 
Run 35 stress 0.1623918 
Run 36 stress 0.1646706 
Run 37 stress 0.1781197 
Run 38 stress 0.1622267 
Run 39 stress 0.1669302 
Run 40 stress 0.1629465 
Run 41 stress 0.1627091 
Run 42 stress 0.1760835 
Run 43 stress 0.1538804 
... Procrustes: rmse 0.01791116  max resid 0.07567211 
Run 44 stress 0.1603052 
Run 45 stress 0.1666163 
Run 46 stress 0.1538248 
... Procrustes: rmse 0.0003646799  max resid 0.001890428 
... Similar to previous best
*** Best solution repeated 1 times

Call:
metaMDS(comm = df1.jac, k = 2, try = 40, trymax = 80, autotransform = FALSE,      maxit = 200) 

global Multidimensional Scaling using monoMDS

Data:     df1.jac 
Distance: binary jaccard 

Dimensions: 2 
Stress:     0.1538244 
Stress type 1, weak ties
Best solution was repeated 1 time in 46 tries
The best solution was from try 1 (random start)
Scaling: centring, PC rotation, halfchange scaling 
Species: scores missing

Run 0 stress 0.09466813 
Run 1 stress 0.09466799 
... New best solution
... Procrustes: rmse 0.0005113888  max resid 0.002601587 
... Similar to previous best
Run 2 stress 0.1116796 
Run 3 stress 0.09466791 
... New best solution
... Procrustes: rmse 5.140713e-05  max resid 0.0001769429 
... Similar to previous best
Run 4 stress 0.09466788 
... New best solution
... Procrustes: rmse 6.431395e-05  max resid 0.0002594411 
... Similar to previous best
Run 5 stress 0.09466814 
... Procrustes: rmse 0.0001741567  max resid 0.0005987574 
... Similar to previous best
Run 6 stress 0.109059 
Run 7 stress 0.09466798 
... Procrustes: rmse 8.15394e-05  max resid 0.0003917043 
... Similar to previous best
Run 8 stress 0.09466792 
... Procrustes: rmse 8.950007e-05  max resid 0.0003379089 
... Similar to previous best
Run 9 stress 0.09466804 
... Procrustes: rmse 0.0001189413  max resid 0.0003712925 
... Similar to previous best
Run 10 stress 0.1004901 
Run 11 stress 0.09466837 
... Procrustes: rmse 0.0002451238  max resid 0.001009793 
... Similar to previous best
Run 12 stress 0.09694634 
Run 13 stress 0.09466805 
... Procrustes: rmse 0.0001311181  max resid 0.0006621851 
... Similar to previous best
Run 14 stress 0.09466838 
... Procrustes: rmse 0.0002467991  max resid 0.00110516 
... Similar to previous best
Run 15 stress 0.1097138 
Run 16 stress 0.09466796 
... Procrustes: rmse 0.0001191739  max resid 0.0005138992 
... Similar to previous best
Run 17 stress 0.09466843 
... Procrustes: rmse 0.0005603326  max resid 0.002919388 
... Similar to previous best
Run 18 stress 0.1005045 
Run 19 stress 0.09466785 
... New best solution
... Procrustes: rmse 4.336385e-05  max resid 0.0001341018 
... Similar to previous best
Run 20 stress 0.09466823 
... Procrustes: rmse 0.0004406253  max resid 0.00239166 
... Similar to previous best
Run 21 stress 0.09466821 
... Procrustes: rmse 0.0002307082  max resid 0.0007713212 
... Similar to previous best
Run 22 stress 0.09693811 
Run 23 stress 0.09466805 
... Procrustes: rmse 0.0001145164  max resid 0.0005228285 
... Similar to previous best
Run 24 stress 0.09466802 
... Procrustes: rmse 0.0001210626  max resid 0.0004425707 
... Similar to previous best
Run 25 stress 0.1005046 
Run 26 stress 0.09466798 
... Procrustes: rmse 0.000102713  max resid 0.0004059971 
... Similar to previous best
Run 27 stress 0.09466796 
... Procrustes: rmse 8.304787e-05  max resid 0.000343865 
... Similar to previous best
Run 28 stress 0.1090607 
Run 29 stress 0.1097141 
Run 30 stress 0.1131481 
Run 31 stress 0.09693847 
Run 32 stress 0.1131483 
Run 33 stress 0.09466825 
... Procrustes: rmse 0.0002105948  max resid 0.0008409997 
... Similar to previous best
Run 34 stress 0.1116801 
Run 35 stress 0.09466791 
... Procrustes: rmse 0.0001849118  max resid 0.0006962707 
... Similar to previous best
Run 36 stress 0.09466793 
... Procrustes: rmse 0.0001164267  max resid 0.0003992448 
... Similar to previous best
Run 37 stress 0.09693871 
Run 38 stress 0.09466835 
... Procrustes: rmse 0.0002585436  max resid 0.001174372 
... Similar to previous best
Run 39 stress 0.1005044 
Run 40 stress 0.09693812 
*** Best solution repeated 11 times

Call:
metaMDS(comm = df1.jac, k = 3, try = 40, trymax = 80, autotransform = FALSE,      maxit = 200) 

global Multidimensional Scaling using monoMDS

Data:     df1.jac 
Distance: binary jaccard 

Dimensions: 3 
Stress:     0.09466785 
Stress type 1, weak ties
Best solution was repeated 11 times in 40 tries
The best solution was from try 19 (random start)
Scaling: centring, PC rotation, halfchange scaling 
Species: scores missing

Again, stress OK with 2D, good with 3D. Look at 2D first

Permutation test for adonis under reduced model
Terms added sequentially (first to last)
Permutation: free
Number of permutations: 999

adonis2(formula = df1.jac ~ type, data = df, permutations = 999)
         Df SumOfSqs      R2      F Pr(>F)
type      2   0.4525 0.04389 1.0328  0.404
Residual 45   9.8587 0.95611              
Total    47  10.3112 1.00000              

No separation seen on Presence-Absence

C

Griffen et al. (2019) examined the oral microbiome of humans, with a focus on the effects of HIV and AntiRetroviral therapy (ART) on this microbiome. They described the microbiome of 341 patients, who fell into one of two categories: HIV^-, and HIV+ with ART. They also matched their samples as far as possible for sex, along with other conditions that can influence oral microbiomes (Candida infection, current smoking). We’ll use their records for these other categories as well. There were more HIV+ than - subjects, but within these groups, approximately equal numbers of two sexes, two Candida infection status, and two smoking categories.

The bacteriome was recorded separately using 16S RNA sequencing, which overed just over 600 taxa.

The data are available from their paper, as two Excel files in the supplementary information. The metadata gives HIV status, and a raft of demographic information. We’ll just work with HIV, sex, Candida, and current smoking. A second sheet has the bacteriome. The code chunk below reads this file (from a Downloads folder - you’ll need to modify the file location). It also screens the file for any bacterial taxa that were not present in this particular study, which reduces the bacterial taxa to 599.

New names:

Use a dissimilarity based approach to assess the combined effects of HIV-ART, sex, Candida infection, and current smoking on the bacteriome

Do these factors act independently on the bacteriome?

Which effects are largest (use some graphical methods)?

First steps: Think about standardization and which distance measure you’ll use. Bray-Curtis seems common for these kind of data, and and counts of different taxa vary widely.

Run factorial model using permanova, based on B-C

Permutation test for adonis under reduced model
Terms added sequentially (first to last)
Permutation: free
Number of permutations: 999

adonis2(formula = df1s.bc ~ HIV * candida * gender * smoking_current, data = df, permutations = 999)
                                    Df SumOfSqs      R2      F Pr(>F)    
HIV                                  1    1.369 0.01250 4.3776  0.001 ***
candida                              1    1.461 0.01333 4.6694  0.001 ***
gender                               1    0.492 0.00449 1.5737  0.016 *  
smoking_current                      1    1.111 0.01014 3.5513  0.001 ***
HIV:candida                          1    0.318 0.00290 1.0151  0.375    
HIV:gender                           1    0.354 0.00324 1.1332  0.214    
candida:gender                       1    0.257 0.00234 0.8202  0.830    
HIV:smoking_current                  1    0.269 0.00245 0.8584  0.773    
candida:smoking_current              1    0.344 0.00314 1.0989  0.250    
gender:smoking_current               1    0.394 0.00360 1.2602  0.096 .  
HIV:candida:gender                   1    0.349 0.00319 1.1164  0.235    
HIV:candida:smoking_current          1    0.280 0.00255 0.8943  0.679    
HIV:gender:smoking_current           1    0.265 0.00242 0.8463  0.784    
candida:gender:smoking_current       1    0.254 0.00232 0.8110  0.873    
HIV:candida:gender:smoking_current   1    0.362 0.00330 1.1562  0.198    
Residual                           325  101.664 0.92809                  
Total                              340  109.541 1.00000                  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Model shows main effects, but no interactions important (or at least, detected)

Could you fit a simpler model to the data?

Run simpler model

Permutation test for adonis under reduced model
Terms added sequentially (first to last)
Permutation: free
Number of permutations: 999

adonis2(formula = df1s.bc ~ HIV + candida + gender + smoking_current, data = df, permutations = 999)
                 Df SumOfSqs      R2      F Pr(>F)    
HIV               1    1.369 0.01250 4.3775  0.001 ***
candida           1    1.461 0.01333 4.6693  0.001 ***
gender            1    0.492 0.00449 1.5736  0.009 ** 
smoking_current   1    1.111 0.01014 3.5512  0.001 ***
Residual        336  105.108 0.95953                  
Total           340  109.541 1.00000                  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

No change; not too surprising, as all four main effects were detected easily.

Now go to data visualization

Run 0 stress 0.1573956 
Run 1 stress 0.1575191 
... Procrustes: rmse 0.008589029  max resid 0.09901144 
Run 2 stress 0.1578788 
... Procrustes: rmse 0.01041235  max resid 0.1078867 
Run 3 stress 0.157925 
Run 4 stress 0.1575481 
... Procrustes: rmse 0.01188138  max resid 0.1520687 
Run 5 stress 0.1579133 
Run 6 stress 0.1576591 
... Procrustes: rmse 0.009942152  max resid 0.1527314 
Run 7 stress 0.157497 
... Procrustes: rmse 0.01111373  max resid 0.1522922 
Run 8 stress 0.1579538 
Run 9 stress 0.1576133 
... Procrustes: rmse 0.01189532  max resid 0.1519556 
Run 10 stress 0.1577343 
... Procrustes: rmse 0.005452067  max resid 0.07804308 
Run 11 stress 0.1578963 
Run 12 stress 0.1579279 
Run 13 stress 0.1576016 
... Procrustes: rmse 0.0109762  max resid 0.107016 
Run 14 stress 0.1579049 
Run 15 stress 0.1589585 
Run 16 stress 0.1578044 
... Procrustes: rmse 0.01067844  max resid 0.1064934 
Run 17 stress 0.1579047 
Run 18 stress 0.158098 
Run 19 stress 0.1572868 
... New best solution
... Procrustes: rmse 0.008393785  max resid 0.1530693 
Run 20 stress 0.1575303 
... Procrustes: rmse 0.007358023  max resid 0.09985515 
Run 21 stress 0.1590169 
Run 22 stress 0.1574836 
... Procrustes: rmse 0.007062616  max resid 0.09973642 
Run 23 stress 0.1589497 
Run 24 stress 0.1575418 
... Procrustes: rmse 0.01043605  max resid 0.1071247 
Run 25 stress 0.1575948 
... Procrustes: rmse 0.007994064  max resid 0.100641 
Run 26 stress 0.1589119 
Run 27 stress 0.157732 
... Procrustes: rmse 0.01101245  max resid 0.1526594 
Run 28 stress 0.157637 
... Procrustes: rmse 0.009312799  max resid 0.1014359 
Run 29 stress 0.1578533 
Run 30 stress 0.1580676 
Run 31 stress 0.1576151 
... Procrustes: rmse 0.009459453  max resid 0.152878 
Run 32 stress 0.1577917 
Run 33 stress 0.157414 
... Procrustes: rmse 0.006382718  max resid 0.09822517 
Run 34 stress 0.1576762 
... Procrustes: rmse 0.01303633  max resid 0.1519725 
Run 35 stress 0.158515 
Run 36 stress 0.1591759 
Run 37 stress 0.1581336 
Run 38 stress 0.1576021 
... Procrustes: rmse 0.004976543  max resid 0.08068629 
Run 39 stress 0.157918 
Run 40 stress 0.1579179 
Run 41 stress 0.1579341 
Run 42 stress 0.1582696 
Run 43 stress 0.157666 
... Procrustes: rmse 0.00894034  max resid 0.1018467 
Run 44 stress 0.157331 
... Procrustes: rmse 0.005409824  max resid 0.09767552 
Run 45 stress 0.1577833 
... Procrustes: rmse 0.006683058  max resid 0.08363389 
Run 46 stress 0.1575939 
... Procrustes: rmse 0.008116432  max resid 0.1084878 
Run 47 stress 0.1575258 
... Procrustes: rmse 0.01301618  max resid 0.1519787 
Run 48 stress 0.1592999 
Run 49 stress 0.1579052 
Run 50 stress 0.1576273 
... Procrustes: rmse 0.01129313  max resid 0.1523282 
Run 51 stress 0.1578542 
Run 52 stress 0.1586214 
Run 53 stress 0.1581151 
Run 54 stress 0.1575617 
... Procrustes: rmse 0.004502256  max resid 0.07990366 
Run 55 stress 0.158704 
Run 56 stress 0.1580142 
Run 57 stress 0.1579508 
Run 58 stress 0.1572867 
... New best solution
... Procrustes: rmse 0.0001807269  max resid 0.002841007 
... Similar to previous best
*** Best solution repeated 1 times

Call:
metaMDS(comm = df1s.bc, k = 3, try = 40, trymax = 80, autotransform = FALSE,      maxit = 200) 

global Multidimensional Scaling using monoMDS

Data:     df1s.bc 
Distance: bray 

Dimensions: 3 
Stress:     0.1572867 
Stress type 1, weak ties
Best solution was repeated 1 time in 58 tries
The best solution was from try 58 (random start)
Scaling: centring, PC rotation, halfchange scaling 
Species: scores missing

Run 0 stress 0.2249552 
Run 1 stress 0.2254015 
... Procrustes: rmse 0.01352051  max resid 0.2271372 
Run 2 stress 0.2255455 
Run 3 stress 0.2258301 
Run 4 stress 0.2253491 
... Procrustes: rmse 0.01645725  max resid 0.2296993 
Run 5 stress 0.22481 
... New best solution
... Procrustes: rmse 0.006257423  max resid 0.08115346 
Run 6 stress 0.226546 
Run 7 stress 0.2270536 
Run 8 stress 0.2255036 
Run 9 stress 0.2253569 
Run 10 stress 0.2248483 
... Procrustes: rmse 0.004470422  max resid 0.08182654 
Run 11 stress 0.2471961 
Run 12 stress 0.2258085 
Run 13 stress 0.2323029 
Run 14 stress 0.2248951 
... Procrustes: rmse 0.004402054  max resid 0.08037907 
Run 15 stress 0.2251572 
... Procrustes: rmse 0.01045667  max resid 0.1510096 
Run 16 stress 0.2253558 
Run 17 stress 0.2254108 
Run 18 stress 0.2250959 
... Procrustes: rmse 0.009892968  max resid 0.1562671 
Run 19 stress 0.225423 
Run 20 stress 0.2248147 
... Procrustes: rmse 0.0005021359  max resid 0.007253184 
... Similar to previous best
Run 21 stress 0.2253845 
Run 22 stress 0.224894 
... Procrustes: rmse 0.004397181  max resid 0.08032982 
Run 23 stress 0.2255606 
Run 24 stress 0.2260425 
Run 25 stress 0.2261645 
Run 26 stress 0.2253466 
Run 27 stress 0.2277116 
Run 28 stress 0.2255894 
Run 29 stress 0.2250862 
... Procrustes: rmse 0.009777865  max resid 0.1540782 
Run 30 stress 0.2257552 
Run 31 stress 0.2253727 
Run 32 stress 0.225797 
Run 33 stress 0.2248478 
... Procrustes: rmse 0.004466226  max resid 0.08183205 
Run 34 stress 0.2250846 
... Procrustes: rmse 0.009723839  max resid 0.1530552 
Run 35 stress 0.2258302 
Run 36 stress 0.2257291 
Run 37 stress 0.22485 
... Procrustes: rmse 0.004479105  max resid 0.08182191 
Run 38 stress 0.2255518 
Run 39 stress 0.2253414 
Run 40 stress 0.2256079 
*** Best solution repeated 1 times

Call:
metaMDS(comm = df1s.bc, k = 2, try = 40, trymax = 80, autotransform = FALSE,      maxit = 200) 

global Multidimensional Scaling using monoMDS

Data:     df1s.bc 
Distance: bray 

Dimensions: 2 
Stress:     0.22481 
Stress type 1, weak ties
Best solution was repeated 1 time in 40 tries
The best solution was from try 5 (random start)
Scaling: centring, PC rotation, halfchange scaling 
Species: scores missing

3D model fits better than 2D. Stress for 2D above 0.2

Start with plots where symbol colour indicates levels of one of the factors

Try some visualizations of pairs of factors using MDS1 & MDS2 (while 3D is preferred, stress of 3D is around .15, vs 0.22, so not huge difference) Separate HIV+/-, then use symbol colour to separate second factor.

References

Dixon, Kelly M., Geoffrey J. Cary, Graeme L. Worboys, and Philip Gibbons. 2018. “The Disproportionate Importance of Long‐unburned Forests and Woodlands for Reptiles.” Ecology and Evolution 8 (22): 10952–63. https://doi.org/gfs587.

Griffen, Ann L., Zachary A. Thompson, Clifford J. Beall, Elizabeth A. Lilly, Carolina Granada, Kelly D. Treas, Kenneth R. DuBois, et al. 2019. “Significant Effect of HIV/HAART on Oral Microbiota Using Multivariate Analysis.” Scientific Reports 9 (1): 19946. https://doi.org/gsfpz3.

Hutto, David, and Kyle Barrett. 2021. “Do Urban Open Spaces Provide Refugia for Frogs in Urban Environments?” Plos One 16 (1). https://doi.org/gr2r9n.